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Computing the lowest eigenvalues of the Fermion matrix by subspace 
iterations 

B. Bunk a * 

a Institut fur Physik, Humboldt-Universitat zu Berlin, Invalidenstr.110, D-10115 Berlin, Germany 

Subspace iterations are used to minimise a generalised Ritz functional of a large, sparse Hermitean matrix. In 
this way, the lowest m eigenvalues are determined. Tests with 1 < m < 32 demonstrate that the computational 
cost (no. of matrix multiplies) does not increase substantially with m. This implies that, as compared to the case 
of a m = 1, the additional eigenvalues are obtained for free. 



1. Introduction 

For a large, sparse Hermitean matrix A G 
(jNxN minimisation of the Ritz functional 



q(x) 



X'X 



xeC 



N 



(1) 



provides the smallest eigenvalue, together with 
the corresponding eigenvector. The minimum of a 
functional can be found iteratively by the method 
of conjugate gradients (CG)|(l[], its application to 
the special case of the Ritz functional has been 
worked out by Geradinj2| and Fried j3j . 

If a small number of lowest eigenvalues is re- 
quired instead, q(x) may be minimised repeat- 
edly, restricting x to the space orthogonal to the 
previous eigenvectors. An efficient implementa- 
tion of this idea is described in p] . 

2. Subspace iterations 

An alternative approach to compute the to 
smallest eigenvalues (with eigenvectors) simulta- 
neously considers the corresponding subspace as 
a whole. An m-dimensional subspace of C N is 
spanned by to non-degenerate (column) vectors, 
which are combined into a rectangular matrix 
x G (j Nxm . From now on, N is considered 'large', 
10 5 say, but m is 'small', e.g. O(10). 

The projector onto the subspace is 



P(x) — x (x'x) 1 x^ . 



(2) 



Note that the inversion of 



G C r ' 



is a 



small problem and can be solved by standard 
techniques. Minimisation is now applied to the 
generalised Ritz functional^ 



q(x) = Tr{P(x)A}. 



(3) 
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At its minimum, P(x) projects onto the space of 
the to lowest eigenvalues. 

An iterative CG minimisation of q(x) starts 
with the computation of the gradient at x: 

g={l-P(x))Ax(x t x)- 1 eC Nxm . (4) 

By construction, x^g — 0. A steepest descent 
method would search a new minimum along g, 
but CG algorithms use an ansatz with a more 
general shift 'direction' h G C Nxm (to be dis- 
cussed later): 

x = x + ha (5) 

with a (small) coefficient matrix a G C mxm . It 
is helpful to assume 

x^h = 0, (6) 

i.e. all shift vectors (columns of h) are orthogonal 
to the current subspace (spanned by the columns 
of x). The matrix a is to be determined to min- 
imise q(x + ha). To bring this to a manageable 
form, consider the projection of A into the space 
spanned by the columns of x and h. This de- 
fines a (2m)-dimensional eigenvalue problem (in 
a non-orthogonal basis), which is to be solved for 
the to lowest eigenvalues. The Jacobi method^] 
can deal with such a 'small' problem. In a matrix 
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notation using m x m blocks, the result can be 
expressed as 



x^ Ax Ah 
h)Ax h)Ah 
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(7) 



The eigenvalues are in A = diag(Ai . . . A m ) and 
the eigenvectors are the columns of 



01 

02 



^ £f 2 m X rn 



(8) 



The new minimising subspace is spanned by 
x0\ + ft./?2- Comparision with eq.(|^) reveals that 
the minimum of q(x + ha) is found at 



a = 020i 1 ■ 

3. Update of the search space 



(9) 



Following the general CG scheme, one starts 
with h — <?, but subsequently h is updated to a 
superposition of the actual gradient and the pre- 
vious search direction. An appropriate ansatz is 



h' = g' + (1 - P{x')) hrt 



(10) 



(g' is the gradient evaluated at x'). The projec- 
tion in front of h insures the orthogonality eq. (^) 
for h! . 7 g (jrnxm j g ano th e r (small) matrix. Its 
choice is somewhat arbitrary: in the linear case 
(minimisation of a quadratic form), it is deter- 
mined by the requirement that the search vectors 
should be conjugate, but this argument does not 
carry over to the nonlinear situation. Experience 
has shown that the following form 



7=(3 t .9)- 1 {5'V-5V} 



(11) 



is reasonably efficient and robust. It is a matrix 
version of the Polak-Ribiere prescription [IjJ. Note 



that matrices like g^g 6 C 



are small. 



4. Convergence and stopping 

By construction, q(x) will decrease monotoni- 
cally with the iteration number n, but the rate of 
convergence and the behavior of individual eigen- 
values is fairly irregular and errors have to be es- 
timated with care. As an example, Fig. [l] shows 



Figure 1. Relative shifts of the four lowest eigen- 
values 

cg28 8 4 random k = .250000 



1 CP 



10" 



l 10 



10" H - 




- 



0.2-10 4 0.4-10 4 0.6-10 4 0.8-10 4 1.0-10 4 
iteration n 



the relative shifts of the four lowest eigenvalues 
of the fermion matrix MtMf for Wilson fermions 
on an 8 4 random SU(2) gauge configuration. 

The following procedure is proposed: assuming 
geometric convergence 

q(x^) » qix^) + af n , 



(12) 



/ is computed from three widely separated n's. 
Then the relation 



q(xW)-q(x(°°>) 



(oo)l 



1 



1-/ 



[q(x ( -^)-q(x ( - n+1 ^](13) 



shows how 'local' shifts are to be rescaled to pro- 
vide error estimates: an eigenvalue A^ is consid- 
ered converged within relative error 5 if 



(14) 



The corresponding eigenvector is frozen and only 
used in further iterations to keep the remaining 
subspace orthogonal to it. 

For the example given above, Fig. [2] demon- 
strates that this procedure leads to correct stop- 
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Figure 2. Relative errors of the four lowest eigen- 
values 
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Table 1 

Number of Av operations to obtain m lowest 
eigenvalues with precision 6 = 10~ 6 on three ran- 




0.2-10 0.4-10 0.6-10 0.8-10 1.0-10 
iteration n 



ping as the eigenvalues approach the precise val- 
ues within a prescribed relative error of S = 1CP 6 . 

5. Computational cost 

Subspace iterations require to store the large 
arrays x, h, Ax, and one more work matrix, this 
amounts to AmN variables in total. 

As with many sparse matrix algorithms, the 
number of multiplications Av (with a single col- 
umn vector v G C ) is considered a fair measure 
for the cost in cpu time. To estimate the per- 
formance of the subspace iterations, tests with 
realistic matrices have been performed: Wilson 
fcrmions in four dimensions with random SU(2) 
gauge fields and (critical) hopping parameter n = 
0.25. In this case, A ~ M^Mf. For lattice sizes 
4 4 ,6 4 ,8 4 , m = 1,2,4,8,16,32 eigenvalues and a 
relative error bound S = 10 -6 , the numbers of op- 
erations Av needed for all eigenvalues to converge 
are shown in Table [l| 

The absolute number of iterations needed to 
achieve a given precision grows with the matrix 
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5268 


21494 


64801 


32 


6000 


25089 


66994 



size (density of the eigenvalues), but the behav- 
ior is difficult to understand in detail. Since di- 
agonalisation within the subspace is done non- 
iteratively, one might speculate that the level 
spacings around the m-th eigenvalue are crucial. 
This requires further studies. 

In any event, the counts in Table [l] show one en- 
couraging feature: the number of Av operations 
does not grow substantially with m. The Av part 
is due to one computation of Ah per iteration 
and costs ~ Nm. It will be dominated for larger 
m by the computation of small matrices like x^x 
(~ Nm 2 per iteration), but these scalar products 
are much cheaper and can be neglected for the 
moderate values of m considered. 
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